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Abstract 

The surface plasmon in simple metal clusters is red-shifted from the Mie 
frequency, the energy shift being significantly larger than the usual spill-out 
correction. Here we develop a variational approach to the RPA collective 
excitations. Using a simple trial form, we obtain analytic expressions for the 
energy shift beyond the spill-out contribution. We find that the additional 
red shift is proportional to the spill-out correction and can have the same 
order of magnitude. 

I. INTRODUCTION 

Simple metal clusters exhibit a strong peak in their optical response that corresponds to a 
collective oscillation of the valence electrons with respect to a neutralizing positively charged 
background. Classically, the frequency of the oscillation is given by the Mie resonance 
formula [1,2], 



where n is the density of a homogeneous electron gas. Quantum finite size effects lead to 
a red shift of this frequency as well as to a redistribution of the oscillator strength (/) into 
closely lying dipole states. Moments of the oscillator strength distribution Mk = Z)j c^f -1 /* 
provide useful information. The first moment Ml, which measures the integral of the f- 
distribution, equals the number of electrons (Thomas- Reiche-Kuhn sum rule). The mean 
square frequency (u 2 ) = M3/M1 is given by the overlap integral of the positive ionic charge 
distribution and the exact ground state electronic density [2]. Within an ionic background 
approximated by a jellium sphere, the mean square frequency is thus exactly related to the 
square Mie frequency by 



where AN/N is the fraction of electrons in the ground state that is outside the jellium sphere 
radius. We called the corresponding energy shift Au so ("spill-out"): 



The actual red shifts are considerably larger than this. For sake of illustrating the dis- 
cussion let us consider the sodium cluster Na^i for which detailed photoabsorption data is 
available [3-7]. The Mie frequency is at 3.5 eV, taking the density corresponding to r s = 3.93 
a.u., while the measured resonance is a peak 2.65 eV having width of about 0.3 eV (FWHM). 
Thus there is a red shift of 24%, which may be compared with a 9% red shift predicted by 
eq. (3) using jellium wave functions. To a large extent clusters with a "magic" number of 
valence electrons behave optically as close shell spherical jellium spheres. The experimen- 
tal photoabsorption spectra for these clusters are well described within the linear response 
theory using either the time-dependent local-density approximation (TDLDA) [8-10] or the 
random phase approximation with exact exchange (RPAE) [11,12]. Red shifts of 14% and 
18% are predicted by time-dependent density functional theory [13] and by the random 
phase approximation [11,12], respectively. The oscillator strength distributions in the RPA 
calculations are typically dominated by a few close states that exhaust almost all of the sum 
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rule. It is this concentration of strength, which we can identify as a dipole surface plasmon, 
that will be of interest in this paper. It is worth of note that singling out a collective state is 
not always possible even in small clusters. Whenever the collective state lies within a region 
of high level density there is a strong fragmentation into p-h states (Landau damping) and 
several excited states may share evenly the strength. We will deal with this problem of the 
definition of the collective state later by proposing a model in which there is no particle-hole 
fragmentation. 

Anharmonic effects in metallic clusters were studied recently by Gerchikov, et ai, [14] 
making use of a coordinate transformation to separate center of mass (cm.) and intrinsic 
motion. The authors show that in absence of coupling between cm. motion and intrinsic 
excitations the surface plasmon associated with a jellium sphere has a single peak which is 
red-shifted with respect to the Mie frequency by the spill-out electrons, Eq. (3). Turning 
on the coupling yields a further red shift which indeed is larger in magnitude than the 
spill-out contribution. Concomitantly, there is a partial transfer of strength into states of 
higher energy preserving the sum rule, Eq. (2). The approach requires the spectrum of 
excitations in the intrinsic coordinates, which were obtained by projection on the computed 
wave functions of the numerical RPAE. 

Another interesting approach to the coupling between the collective and noncollective 
degrees of freedom was developed by Kurasawa, et ai, [15], following the Tomonaga expan- 
sion of the Hamiltonian. The collective coordinate is taken as the cm coordinate, as in ref. 
[14], and the coefficients of the harmonic terms in the Hamiltonian yield Eq. (2) for the 
frequency. The authors derive expressions for the coupling terms in the Hamiltonian and 
use them to estimate the variance of the Hamiltonian in the collective state. They find 
that the variance decreases with size of the cluster as 1/R, where R is the radius of the ion 
distribution. Both the width of the Mie and its shift are obviously related to the variance 
of H, but further assumptions are needed to make a quantitative connection. 

In the present paper we wish to find an analytic estimate of the red shift, keep as far 
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as possible the ordinary formulation of RPA, and not singling out a collective state in the 
Hamiltonian. Our approach will be a variational RPA theory, which we present in the next 
section. The rest of the paper is organized as follows. In Section III we apply the formalism 
to a system of interacting electrons. The model Hamiltonian describes interacting electrons 
confined in a pure harmonic potential, whereas the perturbation corrects for the jellium 
confinement. The model RPA solution is derived analytically and first and second order 
corrections of the frequency shift are given. 



II. VARIATIONAL RPA 

In this section, we establish our notation for the RPA theory of excitations and develop 
a variational expression for perturbations to the collective excitation frequency. The per- 
turbation behaves somewhat differently in RPA than in conventional matrix Hamiltonians 
because the RPA operator is not Hermitean. 

As usual, the starting point is a mean field theory whose ground state is represented by 
an orbital set fa satisfying the orbital equations 

h[po]fa = Cifa (4) 

where p = J2i l<M r )| 2 - The RPA equations are obtained by considering small deviations 
from the ground state, 

< /,.^ < j,. + X(x i e- iut + y i e iut ). (5) 

Here x^yi are vectors in whatever space (r-space, orbital occupation number,...) is used to 
represent fa. The RPA equations can be expressed as 

Sh 

(h[po\ - 6i)xi + 5p* — * fa = uxi (6) 
dp 

-(h[p \ - €i)yi - 5p*^ *fa = ujy { 

bp 
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where the transition density 5p is defined by 

s p = E&te + y») 

i 

and the symbol * denotes an operator or matrix multiplication. Eq. (6) represents linear 
eigenvalue problem for a nonhermitean operator R and the vector \z) = (#i, yi, #2, 2/2, •••)■ 
We will write the equations compactly as 

R\z) = uj\z). 

For a nonhermitean operator, the adjoint vector (z\ is defined as the eigenvector of the 
adjoint equation, (z\R = u(z\. From the symmetry of R it is easy to see that it is given by 
(z\ = (x 1 ,-y 1 ,x 2 ,-y 2 ,...y. 

We now ask how to construct a perturbation theory starting from the zero-order wave 
function \zo) that is the solution of an unperturbed Rq with eigenfrequency u . If we had 
the complete spectrum of i?o, the perturbation series for R = i? + Ai? could be written 
down in the usual way, 

etc. This is in fact what is done in ref. [14]. However, this requires diagonalizing R which 
in general can only be done numerically. 

Instead we shall estimate the energy perturbation using a variational expression for the 
frequency, 

. (z + Xw\R\z + Xw) 
uj = mm — — : — — , (7) 

w (z + \w\z + \w) 

where \w) is a vector to be specified later and A is to be varied to minimize the expression. 
Carrying out the variation and assuming that the perturbation is small, the value of A at 
the minimum is given by 

A _ (zp\R\w) - uj (z \w) 
(w\Rw) — Ui(w\w) 

and the energy shift is 



u « U0 + M^)- ^\- Ul ^ . (9) 

{w\Rw) — ui{w\w) 

Here, lo x = (z \R\z ) = u + (z \AR\z ). 

The next question is how to choose the perturbation \w). With ordinary Hamiltonians, 
one can construct a two-state perturbation theory using the vector obtained by applying 
AR to the unperturbed vector, \w) = AR\z ). However, we will see in the next section that 
this fails completely for the RPA operator. Instead, we will find that an approximation that 
gives qualitatively acceptable results can be made by taking only the re-component of the 
vector defined by applying AR to \z ). 

III. COLLECTIVE LIMIT OF THE SURFACE PLASMON 

We apply the RPA variational perturbation theory derived in the previous section to the 
surface plasmon of small metal clusters. We write the single particle Hamiltonian as 

h = h + AV(r), (10) 
h 2 1 

ho = -— V 2 + -mw V 2 + v * po, (11) 
2m 2 

where v * po is the mean field potential, 

v*p = J v(r,r')p (r')d 3 r f . (12) 

Here v is the electron-electron interaction, which may contain an exchange-correlation con- 
tribution from density functional theory. In this paper, we throughout use the jellium model 
for the ionic background, and also assume that the ion and the electron densities are both 
spherical, ujq and AV(r) are then given by u = Ze 2 /mR 3 and 



AV(r) 



Ze 2 ( 3 r 2 \ Ze 2 



- -- + 



9(r-R), (13) 



r \ 2 2R 2 ) R 

respectively, R being the sharp-cutoff radius for the ion distribution. 

The RPA equations can be solved exactly for the Mie resonance if h is replaced by ho- 
The solution is 
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yj ' 2N \-z<j>) V2iV™j \^/' 

associated with the eigenfrequency uq. Notice that the eigenfrequency ujq is the same as the 
harmonic oscillation frequency in Eq. (11), agreeing with the Kohn's theorem [16-20]. 

To prove that the collective solution (14) satisfies the RPA equation, we use the following 
identity which results from the Hartree-Fock equation, 

(h-e)(A(f)) = [h,A}(f). (15) 

Here A is any one body operator. This yields 

(ho-e)(z(/>) = --d z (/>, (16) 
{ho - e)(d g (j>) = -(mu)%z + (v * d zPo ))(f). (17) 

In the last step, we used the fact that the interaction v is translationally invariant. Notice 
that the transition density is proportional to d z p for the collective solution (14). The 
second term in Eq. (17) is thus exactly canceled by the residual interaction term in the 
RPA equations, proving that the collective ansatz (14) is indeed the eigenfunction of the 
RPA matrix R Q with the eigenvalue u>o. 

The familiar formula relating the red-shift to the electron spill-out probability can be 
recovered from the expectation value of the original RPA matrix, 

— — : — — = lj + Acj. (18) 
{zo\z ) 

However, the wave function z must be taken with the collective ansatz applied to the 
Hamiltonian h. This is different from the z defined in Eq. (14), which was based on the 
Hamiltonian h . In the following, we have no further use for the original z and we will use 
the same name here. Applying the RPA operator Rtoz , we find 

R\zq) = uj \z ) + \u), (19) 

where u is given by 



,„) = { ^ i . (20) 

The expectation value eq. (18) then reduces to 

AiV 

Au = (zq\u) = -^o^r, (21) 

with 

AiV= r Anr 2 drp (r). (22) 
Eq. (21) is just the well-known spill-out formula, Eq.(3), to the first order in AN/N. 



IV. EVALUATION OF THE INTEGRALS 

We now consider the frequency shift in the second order perturbation. Obvious pos- 
sibilities for the perturbation are wo = (y, x) and u, but we find that neither produces a 
significant energy shift. The problem with u is that the x component is tied to the y com- 
ponent in Eq. (20). In fact, the energetics are such the y perturbation is much less than the 
x perturbation. In order to avoid this undesirable feature, as we mentioned in Sec. II, we 
simply take the x component of u for the perturbation. That is, we use 

. . IdAV ( Z( t>\ dAV ( 0\ 

for the \w) in the variational formula (7). With this perturbed wave function, after perform- 
ing the angular integration, we find the three integrals in the formula to be 

. [mu^A-K 3 dAV , 

(zq\u) = -J-TTT-5- / r dr——p (r) 



2N 3 Jr dr 



1 2tt 2 dAVdp 



+J — / r 2 dr=-^, (24) 

V 2Nmuj 3 Jr dr dr ' v ' 

< u|u > = t/« r dr {—) p ° (r) ' (25) 

(z \Ru) = u {z \u) + (u\u) = u (z \u) - \ —^ (u\u). (26) 

y Zl\mujQ 
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In deriving Eq.(26), we have used Eq.(19). We also need to compute (u\R\u) in order to 
estimate the energy shift. Neglecting the residual interaction in the RPA operator R, this 
is expressed as 

dAV dAV 
^| jR | S )^(__0|/,_ e |__0). (27) 

We use Eq.(15) to evaluate the action of the Hamiltonian h onto the u. This yields 
,dAV „ 1 



(h- 



' 2 dAV\ „f^dAV\ „ 
V ^j +2 ( V ^j' V 



(28) 



dz 2m 

Notice that the first term vanishes for the jellium model (13). We thus finally have 

2m 3 j_r ar ar 

In order to get a simple analytic formula for the energy shift, we estimate Eqs. (24), 
(25), (26), and (30) assuming that the density po m the surface region is given by 



Po 



V) ~ Ae- 2K ^ (r > R), (31) 



with k 2 /2m = e, where e is the ionization energy. In order to simplify the algebra, we also 
expand AV and take the first term, 

^^~-3mu)$(r-R). (32) 
dr 

These approximations lead to the following analytic expressions, 

. ._, „ „ 2 fmw^ f R 3 3R 2 9R 3 

\Zq\u) = AnAmaJr, \ < — H H - 

x 1 1 V 2N \Ak 2 4k 3 8k 4 4k 5 

n e ( R 2 R 3 \1 

(«|fl> = nvAm^t (^ + ?g + -L), (34) 



x 4/t 3 4/t 4 4k 5 
4 1 R 2 R 3 



(m u) = U«A^- + (35) 
Note that with the density (31) the spill-out electron number AN is given by 



AiV = Air A (— + + -^—) . (36) 

Retaining only the leading order of 1/kR, we thus have 

(z \u)=mu J l^^-AN, (37) 

(u\u) = 3m 2 ^o (38) 

(u\u) = -3mu -7^1^, (39) 

(u\R\u) = ^muoAN. (40) 

Substituting these expressions into Eq. (9), we finally obtain 

3 (^\ 2 AN 
W = Ul - 16-B- Uo /e (t) (41) 

This is our main result. Note that the perturbation theory breaks down at e = uo/2. In 
realistic situations discussed in the next section, e is always close to u , and the perturbation 
theory should work in principle. 

V. NUMERICAL COMPARISON WITH THE RPA SOLUTIONS 

To assess the reliability of the variational shifts, we have numerically solved the RPA 
equations for the jellium model, using the computer program JellyRpa [13]. A typical 
spectrum is shown in Fig. 1. This represents Na 2 o as a system of 20 electrons in a background 
spherical charge distribution with a density corresponding to r s = 3.93 a.u. and total charge 
Q = 20. The strength function includes an artificial width of T = 0.1 eV for display purposes. 
The Mie frequency, Eq. (1), is indicated by uj , while the prediction of the spill-out formula, 
Eq.(3), is shown as oo so in the figure. One sees that the strength function is fragmented 
into two large components that are considerably red-shifted from the Mie frequency, and 
smaller contributions at higher frequencies. The corresponding spectrum with the jellium 
background potential replaced by a pure harmonic potential is shown by the dashed line. 
The numerical RPA frequency agrees very well with the Mie value in this case, showing that 
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the numerical algorithms used in JellyRpa are sufficiently accurate for our purposes. The 
red shift can be more easily displayed by a plot of the integrated strength function, shown 
in the lower panel of the figure. If we define the shift as the point where the integrated 
strength reaches half of the maximum value, it corresponds to 5uj = 0.166 c^Mie- On the 
other hand, the collective formula for the red shift, Eq. (3), only gives 5u = 0.058 cj Mic , 
when the integral for AN is evaluated with the ground state density. 

The strength becomes increasingly fragmented in heavier clusters, making a precise def- 
inition of the red shift problematic. We therefore have simplified the jellium model in our 
numerical computations to see the effects of the shift without the fragmentation of the 
strength that occurs physically. To this end we put all the electrons in the lowest s-orbital, 
treating them as bosons. Otherwise, the model is the same as the usual jellium model, with 
the electron orbitals determined self-consistently in a background charge density of a uni- 
form sphere. This model is easily implemented with JellyRpa by assigning the occupation 
probabilities of the orbitals appropriately. Taking the density parameter as r s = 3.93 a.u., 
appropriate for sodium clusters, one finds that the ionization potential is rather close to 
the value of the usual (fermionic) jellium model. For example, in the cluster with N = 20 
atoms, the ionization potential e has a value 2.84 eV for usual jellium model and the value 
4.11 eV for our simplified s-wave treatment. 

The results of the numerical calculation with the full effect of the surface are shown in 
Fig. 2 as the solid line. The collective spill-out correction from Eq. (3) is also shown as 
the dotted line. One sees that the additional shift due to the wave function perturbation 
is comparable to the spill-out correction, and has a similar iV-dependence. The shift given 
by the variational formula Eq. (9) is shown by the dashed line. The functional dependence 
predicted by the formula is confirmed by the numerical calculations, but the coefficient of 
N is too small by a factor of two or so. 
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VI. CONCLUDING REMARKS 



We have developed a variational approach to treat perturbations to the collective RPA 
wave functions, and have applied it to the surface plasmon in small metal clusters. Our 
zeroth order solution is the same as that used by Gerchikov et al. [14] and Kurasawa et 
al. [15]. It corresponds to the center of mass motion, and is the exact RPA solution when 
the ionic background potential is a harmonic oscillator. The deviation of the background 
potential from the harmonic shape is responsible for the perturbation. The first order 
perturbation yields the well-known spill-out formula for the plasmon frequency, as was also 
shown in Refs. [14,15]. The higher order corrections lead to the additional energy shift 
of the frequency [14], the anharmonicity of the spectrum [14], and the fragmentation of 
the strength [15]. Those effects were studied in Refs. [14,15] by considering explicitly the 
couplings between the center of mass and the intrinsic motions. In this paper, we assumed 
some analytic form for the perturbation and determined its coefficient variationally. We 
found that this approach qualitatively accounts for the red shift of the collective frequency, 
but its magnitude came out too small by about a factor of two. 

In order to have a more quantitative result, one would have to improve the variational 
wave function. An obvious way is to introduce more than one term. Our method may 
be viewed as the first iteration of any iterative method for RPA [21-23]. One may need 
more than one iteration to get a convergence and thus a sufficiently large energy shift. 
Another possible way is to construct the perturbed wave function based on the local RPA. 
The authors of Ref. [24] expanded the collective operator with local functions and solved 
a secular equation to determine the frequency. They showed that the expansion of the 
collective operator with three functions, r cos#,r 3 cos6>, and r 5 cos#, gives a satisfactory 
result for the collective frequency. 

The method developed in this paper is general, and is not restricted to the surface 
plasmon in micro clusters. One interesting application may be to the giant dipole resonance 
in atomic nuclei. In heavy nuclei, the mass dependence of the isovector dipole frequency 
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deviates from the prediction of the Goldhaber- Teller model, that is based on a simple cm. 
motion [25,26]. The shift of collective frequency can be attributed to the effect of deviation 
of the mean-field potential from the harmonic oscillator, and a similar treatment as the 
present one is possible. 
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FIG. 1. Strength function of Na2o i n the jellium model. Upper panel shows the dipole strength 
function, broadened by a artifical width. Lower panel shows the integerated strength function. 
Dashed line is the results of the computation in which the jellium background potential is replaced 
by a harmonic oscillator. 
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FIG. 2. Collective excitation frequency in the s-wave jellium model as a function of N. The 
solid line is the result of the numerical calculation. This is compared with the spill-out formula eq. 
(3) and the perturbation formula eq. (9) as the dotted and dashed lines, respectively. 
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